↜ Back to index Introduction to Numerical Analysis 2

Lecture 6

Gaussian elimination (ガウスの消去法)

Our goal is to now solve a linear system Ax = b for any invertible (regular, 正則) matrix A.

We know how to do this if A is upper- or lower triangular (Lecture a4).

We therefore need to convert the system Ax=b into a form with a triangular matrix. We can do this rather easily by performing row operations on A and b. Indeed, if we add a multiple of an one equation to another equation in the system, we do not change the solution. This is equivalent to adding a multiple of a row of the matrix A to another row, and adding the multiple of the component of b in the same row to the component of b in the other row.

Example 1. For the following A and b, \tag{1} \begin{pmatrix} 1 & 3 & 1 \\ 1 & 1 & -1\\ 3 & 11 & 6 \end{pmatrix}, \begin{pmatrix} 9 \\ 1 \\ 34 \end{pmatrix}, we can add (-1) \times {\rm first\ row} to the second row, and (-3) \times {\rm first\ row} to the third row and we get \tag{2} \begin{pmatrix} 1 & 3 & 1 \\ 0 & -2 & -2\\ 0 & 2 & 3 \end{pmatrix}, \begin{pmatrix} 9 \\ -8 \\ 7 \end{pmatrix} We can then add the second row to the third row and obtain \tag{3} \begin{pmatrix} 1 & 3 & 1 \\ 0 & -2 & -2\\ 0 & 0 & 1 \end{pmatrix}, \begin{pmatrix} 9 \\ -8 \\ -1 \end{pmatrix} System with this matrix and right hand side is now easy to solve.

Exercise 1. Write a code that reads n, a and b (the first row of a is read first, then the second row, etc.), and then adds a multiple of the first row of a and b to the other rows so that the all the elements in the first column but the top one are 0.

The program should print the resulting array a as a matrix and the resulting b.

To print the array a as a matrix, you can use the code from the previous lecture:
```
do i = 1, m
    print *, a(i, :)
end do
```

As a test, download the file data.txt, save it in the directory where you run the program, and then run your program as

./a.exe < data.txt

data.txt contains A and b as in (1) above. The output therefore should be as in (2) above:

$ ./a.exe < data.txt                     
 a =
   1.00000000       3.00000000       1.00000000    
   0.00000000      -2.00000000      -2.00000000    
   0.00000000       2.00000000       3.00000000    
 b =
   9.00000000    
  -8.00000000    
   7.00000000

Exercise 2. Modify the program in exercise 1 so that it performs the operation on all columns except the last one.

That is, for column 2, add a multiple of the second row to rows 3, …, n so that all the elements in column 2 on rows 3, …, n are 0. Continue for columns 3, …, n-1.

Using data.txt, you should get an output as in (3) above.

$ ./a.exe < data.txt 
 a =
   1.00000000       3.00000000       1.00000000    
   0.00000000      -2.00000000      -2.00000000    
   0.00000000       0.00000000       1.00000000    
 b =
   9.00000000    
  -8.00000000    
  -1.00000000

Download data2.txt which contains a 5\times 5 matrix. The output should be

$ ./a.exe < data2.txt           
 a =
   1.00000000       5.00000000       4.00000000       6.00000000       2.00000000    
   0.00000000      -4.00000000      -3.00000000      -11.0000000       9.00000000    
   0.00000000       0.00000000      -1.00000000       8.00000000      -20.0000000    
   0.00000000       0.00000000       0.00000000       32.2500000      -63.7500000    
   0.00000000       0.00000000       0.00000000       0.00000000       113.744186    
 b =
   1.00000000    
   1.00000000    
  -1.00000000    
  -3.00000000    
   6.02325535

Exercise 3. Combine the code from Exercise 2 with the code for solving the system with triangular matrix to find the solution x of the linear system Ax = b given a matrix A in an array a and vector b in an array b.

The program should output the following:

$ ./a.exe < data.txt
  -5.00000000    
   5.00000000    
  -1.00000000

$ ./a.exe < data2.txt
   1.63013697    
 -0.188509524    
   3.41444016E-02
   1.16540482E-02
   5.29544018E-02

Partial pivoting

We can solve almost any linear system with an invertible matrix. However, our algorithm fails whenever a(i, i) is 0 when we want to make elements in column i in rows i+1, …, n zero.

Example 2. System with A and b like \begin{pmatrix} 0 & 1\\ 1 & 0 \end{pmatrix}, \begin{pmatrix} 1\\ 2 \end{pmatrix} clearly has a unique solution, but our algorithm in Exercise 3 fails.

Try this on your program from Exercise 3 with data3.txt.

One idea to overcome this issue is to reorder rows when we reach a situation like this.

In partial pivoting, when we want to make entries of the matrix in column k zero in rows k+ 1, …, n, we first perform the following 2 steps:

We find the row m with k \leq m \leq n such that a(m, k) has the largest absolute value.
If k \neq m, we swap the rows k and m of the array a and vector b.

Then we proceed as usual.

Example 3. Consider the system in Example 1: \begin{pmatrix} 1 & 3 & 1 \\ 1 & 1 & -1\\ 3 & 11 & 6 \end{pmatrix}, \begin{pmatrix} 9 \\ 1 \\ 34 \end{pmatrix},

We consider column k = 1.

The row with the largest absolute value of an entry in column 1 is row m = 3.
We swap row 1 and 3 and end up with

\begin{pmatrix} 3 & 11 & 6\\ 1 & 1 & -1\\ 1 & 3 & 1 \end{pmatrix}, \begin{pmatrix} 34\\ 1 \\ 9 \end{pmatrix},

The we proceed in adding (-1/3) \times {\rm first\ row} to both second and third rows.

Exercise 4. Write a code that given:

n a positive integer;
a n\times n matrix A stored in a two-dimensional array a;
a vector b of size n stored in a one-dimensional array b

prints out the solution x_1, x_2, …, x_n of the linear system A x = b

using Gaussian elimination with partial pivoting.

The program should read the input data as in Exercise 1.

The program must print out only:

If when considering column k, none of the entries a(k, k), a(k+ 1, k), …, a(n, k) are nonzero:
```
Error: no nozero entry when considering column k   
```
where k is the actual number of the column.
Otherwise: It prints the values x_1, x_2, …, x_n separated by new lines.

Download the test files:

Expected outputs with the individual test files:

$ ./a.exe < data.txt                        
  -5.00000000    
   4.99999952    
 -0.999999523    
$ ./a.exe < data2.txt
   1.63013697    
 -0.188509509    
   3.41443419E-02
   1.16540594E-02
   5.29544093E-02
$ ./a.exe < data3.txt
   2.00000000    
   1.00000000    
$ ./a.exe < data4.txt
 Error: no nozero entry when considering column            3

Hints.

To find the row with the element with largest absolute value, you can use the algorithm we covered in Lecture a3, Exercise 5 (solution).
To stop a Fortran program early, you can use command stop.
Make sure that there are no “index out of bounds” errors by compiling with -fcheck=bounds option as explained in Lecture a3.
Make sure that the program produces the correct output for the example inputs. More or fewer spaces or new lines are OK, but no extra text in the output!